library(sp)
library(raster)
library(R.matlab)
library(plyr)
library(data.table)
library(maptools)
library(rgdal)
library(spatstat)
setwd("~/Desktop/Fall_2018_Classes/GIS_713/Final_Project")
library(plyr)
library(data.table)
CHIPS_df = read.table("chips.csv", header = TRUE, row.names=NULL, sep=",")
CHIPS_df <- CHIPS_df[!CHIPS_df$lat < 500000.00, ]
CHIPS_df <- CHIPS_df[!CHIPS_df$long < 50000.00, ]
# correcting odd lat coords
CHIPS_df$lat <- CHIPS_df$lat / 10000
CHIPS_df$long <- CHIPS_df$long / 10000
# And a function that shifts vectors conviniently (we gonna need this later):
shift.vec <- function (vec, shift) {
if(length(vec) <= abs(shift)) {
rep(NA ,length(vec))
}else{
if (shift >= 0) {
c(rep(NA, shift), vec[1:(length(vec)-shift)]) }
else {
c(vec[(abs(shift)+1):length(vec)], rep(NA, abs(shift))) } }
}
# Now, let’s calculate the distances between successive positions and the respective speed in this segment.
# Shift vectors for lat and lon so that each row also contains the next position.
CHIPS_df$lat.p1 <- shift.vec(CHIPS_df$lat, -1)
CHIPS_df$lon.p1 <- shift.vec(CHIPS_df$long, -1)
# Calculate distances (in metres) using the function pointDistance from the ‘raster’ package.
# Parameter ‘lonlat’ has to be TRUE!
CHIPS_df$dist.to.prev <- apply(CHIPS_df, 1, FUN = function (row) {
pointDistance(c(as.numeric(row["lat.p1"]), as.numeric(row["long.p1"])),
c(as.numeric(row["lat"]), as.numeric(row["long"])), lonlat = T)
})
# Transform the column ‘time’ so that R knows how to interpret it.
CHIPS_df$time_new <- strptime(CHIPS_df$initial_time_stamp_mat, format="%m/%d/%Y %H:%M")
# Shift the time vector, too.
CHIPS_df$time.p1 <- shift.vec(CHIPS_df$time_new, -1)
# Calculate the number of seconds between two positions.
CHIPS_df$time.diff.to.prev <- as.numeric(difftime(CHIPS_df$time.p1, CHIPS_df$time_new))
# Calculate metres per seconds, kilometres per hour and two LOWESS smoothers to get rid of some noise.
CHIPS_df$speed.m.per.sec <- CHIPS_df$dist.to.prev / CHIPS_df$time.diff.to.prev
CHIPS_df$speed.km.per.h <- CHIPS_df$speed.m.per.sec * 3.6
CHIPS_df$speed.km.per.h <- ifelse(is.na(CHIPS_df$speed.km.per.h), 0, CHIPS_df$speed.km.per.h)
CHIPS_df$lowess.speed <- lowess(CHIPS_df$speed2, f = 0.2)$y
CHIPS_df$lowess.alt <- lowess(CHIPS_df$altitude, f = 0.2)$y
CHIPS_df$lowess.conduct <- lowess(CHIPS_df$conductance_z, f = 0.2)$y
setnames(CHIPS_df, "long", "lon")
pt1 <- CHIPS_df[CHIPS_df$participant == 1, ]
pt2 <- CHIPS_df[CHIPS_df$participant == 2, ]
pt3 <- CHIPS_df[CHIPS_df$participant == 3, ]
pt4 <- CHIPS_df[CHIPS_df$participant == 4, ]
pt5 <- CHIPS_df[CHIPS_df$participant == 5, ]
pt6 <- CHIPS_df[CHIPS_df$participant == 6, ]
pt7 <- CHIPS_df[CHIPS_df$participant == 7, ]
pt8 <- CHIPS_df[CHIPS_df$participant == 8, ]
pt9 <- CHIPS_df[CHIPS_df$participant == 9, ]
pt10 <- CHIPS_df[CHIPS_df$participant == 10, ]
pt11 <- CHIPS_df[CHIPS_df$participant == 11, ]
# Now, let’s plot all the stuff!
# Plot elevations and smoother
layout(matrix(1:3, nrow=3))
plot(CHIPS_df$altitude, type = "l", bty = "n", xaxt = "n", lwd= 3,
ylab = "Elevation",
xlab = "", col = "grey60")
lines(CHIPS_df$lowess.alt, col = "green", lwd = 3)
abline(h = mean(CHIPS_df$altitude), lty = 2, lwd = 3, col = "green")
legend(x="bottomright", legend = c("GPS elevation", "LOWESS elevation",
"Mean elevation"),
col = c("grey60", "green", "green"), lwd = c(2,4,2), lty = c(1,2,2),
bty = "n")
# Plot speeds and smoother
plot(CHIPS_df$speed2, type = "l", bty = "n", lwd= 3, xaxt = "n",
ylab = "Speed (km/h)", xlab = "", col = "grey60")
lines(CHIPS_df$lowess.speed, col = "red", lwd = 3)
abline(h = mean(CHIPS_df$speed2), lty = 2, lwd = 3, col = "red")
legend(x="topright", legend = c("GPS speed", "LOWESS speed",
"Mean speed"),
col = c("grey60", "red", "red"), lwd = c(2,4,2), lty = c(1,2,2),
bty = "n")
# Plot conductnace and smoother
plot(CHIPS_df$conductance_z, type = "l", bty = "n", lwd= 3, xaxt = "n",
ylab = "Skin Conductance", xlab = "", col = "grey60")
lines(CHIPS_df$lowess.conduct, col = "blue", lwd = 3)
abline(h = mean(CHIPS_df$conductance_z), lty = 2, lwd = 3, col = "blue")
legend(x="topright",
legend = c("Conductance", "LOWESS conductance", "Mean conductance"),
col = c("grey60", "blue", "blue"), lwd = c(2,4,2), lty = c(1,2,2),
bty = "n")
par(mar=c(5, 4, 4, 2) + 0.1)
# Plotting the elevation and timestamp of each waypoint using ggplot, allows us to visualise the hike towards the valley of Grauson. As a preparatory step, we use the ymd_hms() function from the lubridate library to convert the string representating the timestamp to a proper R time-object. As to not confuse ggplot, we also do not pass the SpatialPointsDataFrame-object directly, but convert it to a regular dataframe with as.data.frame():
if(!require(lubridate)) install_github("rstudio/lubridate")
if(!require(ggplot2)) install_github("rstudio/ggplot2")
if(!require(gridExtra)) install_github("rstudio/gridExtra")
# plot of time and elevation, colored by skin conductance
time_ele_conduct_plot <- ggplot(as.data.frame(CHIPS_df), # convert to regular dataframe
aes(x=time, y=altitude, color = conductance_z)) +
scale_color_gradient(low="blue", high="red") +
geom_point(alpha = 0.01, size = 4) +
labs(x='Cycling time', y='Elevations (meters)')
# plot of time and speed, colored by skin conductance
time_speed_conduct_plot <- ggplot(as.data.frame(CHIPS_df), # convert to regular dataframe
aes(x=time, y=speed2, color = conductance_z)) +
scale_color_gradient(low="blue", high="red") +
geom_point(alpha = 0.08, size = 4) +
labs(x='Cycling time', y='Speed (km/h)')
grid.arrange(time_ele_conduct_plot, time_speed_conduct_plot, nrow=2)
shp_dsn <- "~/Desktop/Fall_2018_Classes/GIS_713/Final_Project/NL006L3_TILBURG/Shapefiles"
landcover <- readOGR(path.expand(shp_dsn), 'NL006L3_TILBURG_UA2012')
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/garrettmillar/Desktop/Fall_2018_Classes/GIS_713/Final_Project/NL006L3_TILBURG/Shapefiles", layer: "NL006L3_TILBURG_UA2012"
## with 12086 features
## It has 11 fields
# Projection
landcover <- spTransform(landcover, CRS('+proj=longlat'))
# Conversion into SpatialPoints
coordinates(CHIPS_df) <- ~ lon + lat
coordinates(pt1) <- ~ lon + lat
coordinates(pt2) <- ~ lon + lat
coordinates(pt3) <- ~ lon + lat
coordinates(pt4) <- ~ lon + lat
coordinates(pt5) <- ~ lon + lat
coordinates(pt6) <- ~ lon + lat
coordinates(pt7) <- ~ lon + lat
coordinates(pt8) <- ~ lon + lat
coordinates(pt9) <- ~ lon + lat
coordinates(pt10) <- ~ lon + lat
coordinates(pt11) <- ~ lon + lat
# Setting default projection
proj4string(CHIPS_df) <- CRS('+proj=longlat')
proj4string(pt1) <- CRS('+proj=longlat')
proj4string(pt2) <- CRS('+proj=longlat')
proj4string(pt3) <- CRS('+proj=longlat')
proj4string(pt4) <- CRS('+proj=longlat')
proj4string(pt5) <- CRS('+proj=longlat')
proj4string(pt6) <- CRS('+proj=longlat')
proj4string(pt7) <- CRS('+proj=longlat')
proj4string(pt8) <- CRS('+proj=longlat')
proj4string(pt9) <- CRS('+proj=longlat')
proj4string(pt10) <- CRS('+proj=longlat')
proj4string(pt11) <- CRS('+proj=longlat')
library(maptools)
library(rgdal)
library(leaflet)
require(pals)
conduct.pal <- colorNumeric (c("dodgerblue4", "slategray2", "red3"), pt1$conductance_z)
m <- leaflet() %>%
# Add tiles
addProviderTiles("Esri.WorldTopoMap", group = "Topographical") %>%
addProviderTiles("OpenStreetMap.Mapnik", group = "Road map") %>%
addProviderTiles("Esri.WorldImagery", group = "Satellite") %>%
addLegend(position = 'bottomright', opacity = 0.4,
colors = "blue",
labels = 'Participant 1',
title = 'Cycling Routes (Netherlands)') %>%
# Layers control
addLayersControl(position = 'topright',
baseGroups = c("Topographical", "Road map", "Satellite"),
options = layersControlOptions(collapsed = FALSE)) %>%
addCircles (data=pt1, group='Cycling routes', stroke = T, radius = 100,
fillOpacity = 0.2, fillColor = conduct.pal(pt1$conductance_z),
opacity = 0.2, color = conduct.pal(pt1$conductance_z)) %>%
addLegend(values = pt1$conductance_z, pal = conduct.pal,
opacity = 1, title = "Skin Conductivity", position = "bottomleft")
m
urban <- c("Discontinuous dense urban fabric (S.L. : 50% - 80%)",
"Industrial, commercial, public, military and private units",
"Discontinuous low density urban fabric (S.L. : 10% - 30%)",
"Isolated structures",
"Continuous urban fabric (S.L. : > 80%)",
"Discontinuous very low density urban fabric (S.L. : < 10%)",
"Discontinuous medium density urban fabric (S.L. : 30% - 50%)",
"Construction sites",
"Mineral extraction and dump sites")
natural <- c("Arable land (annual crops)",
"Pastures",
"Forests",
"Land without current use",
"Green urban areas",
"Permanent crops (vineyards, fruit trees, olive groves)",
"Sports and leisure facilities",
"Herbaceous vegetation associations (natural grassland, moors...)",
"Open spaces with little or no vegetation (beaches, dunes, bare rocks, glaciers)",
"Wetlands")
roads <- c("Other roads and associated land",
"Railways and associated land",
"Fast transit roads and associated land")
name = landcover$ITEM2012
landcover$group <- with(landcover, ifelse(name %in% urban, "urban",
ifelse(name %in% natural,
"natural", "roads")))
col = terrain.colors(3)
plot(landcover, col = col, col.regions = landcover$group,
edge.col = "transparent", axes = F,
colorkey = list(space = "bottom", height = 0.5,
width = 0.7),
main = "Study Area", sub = "Types of Landcover",
par.settings = list(axis.line = list(col = 'transparent')))
legend("bottomleft", title = "Landcover", fill = col,
legend = c("Urban", "Natural", "Roads"))
# clipping landcover polygon to cycling route
neth_clipped <- as(crop(landcover, CHIPS_df), "SpatialPolygonsDataFrame")
col = terrain.colors(3)
conduct.pal <- colorNumeric (c("dodgerblue4", "slategray2", "red3"), pt1$conductance_z)
layout(matrix(1:6, nrow=1))
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
edge.col = "transparent", axes = F,
colorkey = list(space = "bottom", height = 0.5, width = 0.7),
main = "",
par.settings = list(axis.line = list(col = 'transparent')))
points(pt1, col = conduct.pal(pt1$conductance_z), cex = 1.5, pch = 20 )
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
edge.col = "transparent", axes = F,
colorkey = list(space = "bottom", height = 0.5, width = 0.7),
main = "",
par.settings = list(axis.line = list(col = 'transparent')))
points(pt2, col = conduct.pal(pt2$conductance_z), cex = 1.5, pch = 20 )
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
edge.col = "transparent", axes = F,
colorkey = list(space = "bottom", height = 0.5, width = 0.7),
main = "",
par.settings = list(axis.line = list(col = 'transparent')))
points(pt3, col = conduct.pal(pt3$conductance_z), cex = 1.5, pch = 20 )
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
edge.col = "transparent", axes = F,
colorkey = list(space = "bottom", height = 0.5, width = 0.7),
main = "",
par.settings = list(axis.line = list(col = 'transparent')))
points(pt4, col = conduct.pal(pt4$conductance_z), cex = 1.5, pch = 20 )
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
edge.col = "transparent", axes = F,
colorkey = list(space = "bottom", height = 0.5, width = 0.7),
main = "",
par.settings = list(axis.line = list(col = 'transparent')))
points(pt5, col = conduct.pal(pt5$conductance_z), cex = 1.5, pch = 20)
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
edge.col = "transparent", axes = F,
colorkey = list(space = "bottom", height = 0.5, width = 0.7),
main = "",
par.settings = list(axis.line = list(col = 'transparent')))
points(pt6, col = conduct.pal(pt6$conductance_z), cex = 1.5, pch = 20)
library(OpenStreetMap)
map <- openmap(as.numeric(c(max(pt1$lat), min(pt1$lon))),
as.numeric(c(min(pt1$lat), max(pt1$lon))), type = "stamen-terrain")
transmap <- openproj(map, projection = proj4string(pt1))
map3d <- function(map, ...){
if(length(map$tiles)!=1){stop("multiple tiles not implemented") }
nx = map$tiles[[1]]$xres
ny = map$tiles[[1]]$yres
xmin = map$tiles[[1]]$bbox$p1[1]
xmax = map$tiles[[1]]$bbox$p2[1]
ymin = map$tiles[[1]]$bbox$p1[2]
ymax = map$tiles[[1]]$bbox$p2[2]
xc = seq(xmin,xmax,len=ny)
yc = seq(ymin,ymax,len=nx)
colours = matrix(map$tiles[[1]]$colorData,ny,nx)
m = matrix(0,ny,nx)
surface3d(xc,yc,m,col=colours, ...)
}
library(RColorBrewer)
# display.brewer.all()
bp = brewer.pal(11,"RdBu")
library(colourschemes)
cs = rampInterpolate(pt1$conductance_z, rev(bp))
pt1$conduct_pal <- cs(pt1$conductance_z)
plot3d(pt1$lon, pt1$lat, pt1$time, xlab="Longitude",
ylab="Latitude", zlab="Time", type = "s",
col = pt1$conduct_pal, size = 2.5, nam)
You must enable Javascript to view this page properly.